
disp('-')
disp('-')
disp('-')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%                    TABLE 6:               %%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('-')
disp('-')
disp('-')


pretemp    = pre ;
pre        = pre*36 ;
Postscore  = kappapostEB*36 ; 
DEPVAR_MeanSD = [mean(kappapost*36,'omitnan') std(kappapostEB*36,'omitnan')] %#ok 
Deltakappa_MeanSD = [mean((kappapost-kappapre)*36,'omitnan') std((kappapostEB-kappapreEB)*36,'omitnan')] %#ok 


[Ftemp,Xtemp]    = ecdf(AchievementData.Q(AchievementData.Q>=2)); Ftemp(2)=[]; Xtemp(2)=[] ;
tempbenchmark    = interp1(Ftemp,Xtemp,[(0.5-0.341); 0.5],'linear') ;
BenchmarkChangeQ = tempbenchmark(2)-tempbenchmark(1) %#ok
[Ftemp,Xtemp]    = ecdf(AchievementData.T(AchievementData.Q>=2)) ; Ftemp(2)=[]; Xtemp(2)=[] ; 
tempbenchmark    = interp1(Ftemp,Xtemp,[(0.5-0.341); 0.5],'linear') ;
BenchmarkChangeT = tempbenchmark(2)-tempbenchmark(1) %#ok


Tmean   = mean(T) ;
Tsig    = std(T) ;
T       = (T)/Tsig ; %%%%Here we normalize T for numerical stability
BenchmarkChangeTstd = (BenchmarkChangeT)/Tsig ;

Qmean   = mean(Q) ;
Qsig    = std(Q) ;
Q       = (Q)/Qsig ;
BenchmarkChangeQstd = (BenchmarkChangeQ)/Qsig ;
Premean = mean(pre) ;
Presig  = std(pre) ;
prestdz = (pre-Premean)/Presig ; %%%%We use standardized pre-test score in this section for numerical stability in Model 6
disp('******************************************************************************************************') 
disp('***********   HETEROSKEDASTICITY ROBUST ANALYSIS: SHORT-RUN HC PRODUCTION   **************************') 
disp('******************************************************************************************************') 
StDevDeltaScore  = std((kappapostEB-kappapreEB)*36,'omitnan') %#ok ; % std(Deltascore,'omitnan') ; 
PostscoreMultiplier = StDevDeltaScore/std(Postscore,'omitnan') %#ok

ExpVars_SR1 = [ones(size(Postscore)), ... 1
               T, ... 2
               T.^2, ... 3
               Q, ... 4
               Q.^2, ...5
               T.*Q] ; %6
labels_SR1  = {'Const*';
               'T';'T^2'; ...
               'Q';'Q^2';'T*Q'} ; 
[alpha_SR1,se_SR1,tstat_SR1,pval_SR1,CI_SR1,out_SR1] = regressHet(Postscore,ExpVars_SR1,'FGLSrobust',labels_SR1,[],[],0,[],[],[]); 

StDevEffectATE    = nan(length(alpha_SR1)+2,1) ;
%%%%NOTE: in the following line, recall that std(T)=1 after the normalization.
dDeltadT          = (alpha_SR1(2)) + 2*(alpha_SR1(3)).*T + Q.*(alpha_SR1(6)) ;
StDevEffectATE(2) = mean( BenchmarkChangeTstd*dDeltadT/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD effect of T
dDeltadQ          = (alpha_SR1(4)) + 2*(alpha_SR1(5)).*Q + T.*(alpha_SR1(6)) ;
StDevEffectATE(4) = mean( BenchmarkChangeQstd*dDeltadQ/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD effect of Q
OtherEffect    = nan(length(alpha_SR1)+2,1) ;
%%%%This is a mean total derivative of kappa with respect to Q, for a
%%%%change of BenchmarkChangeQstd, as a fraction of StDevDeltaScore
OtherEffect(4) = mean( ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd/StDevDeltaScore ) ;  %%%%This is a mean total derivative effect of Q (in original DeltaS units)
disp('-----------------------------------------------------------------------------------')
disp('-------------------------  SHORT-RUN PROGRESS: MODEL 1  ---------------------------')
disp('-----------------------------------------------------------------------------------')
baselinetemp = sort( alpha_SR1(1) - pre ) ;
BASELINE_SR1 = [mean(baselinetemp) median(baselinetemp) std(baselinetemp)] %#ok
% BASELINE_SR1 = [alpha_SR1(1) alpha_SR1(1) 0 alpha_SR1(1) alpha_SR1(1)] - kappapreEB*36 %#ok
out_SR1.results.StDevEffectATE = StDevEffectATE;
out_SR1.results.OtherEffect = OtherEffect;
WhiteTest  = out_SR1.White %#ok
Method     = out_SR1.method %#ok
disp(out_SR1.results(:,:))
%%%%Test for joint significance of all T terms:           
R = zeros(3,length(ExpVars_SR1(1,:))) ; R(1,2)=1; R(2,3)=1; R(3,6)=1; r = zeros(3,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR1,'FGLSrobust',labels_SR1,R,r,0,[],[],[]); 
TTest = temp.Wald %#ok
%%%%Test for joint significance of all Q terms:           
R = zeros(3,length(ExpVars_SR1(1,:))) ; R(1,4)=1; R(2,5)=1; R(3,6)=1; r = zeros(3,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR1,'FGLSrobust',labels_SR1,R,r,0,[],[],[]); 
QTest = temp.Wald %#ok    

%%


ExpVars_SR15 = [ones(size(Postscore)), ... 1
               T, ... 2
               T.^2, ... 3
               Q, ... 4
               Q.^2, ...5
               T.*Q, ...6
               prestdz, ...7 %%%%We use standardized pre-test score for numerical stability
               D2, ...8
               D3, ...9
               D2.*T, ...10
               D3.*T, ...11
               D2.*Q, ...12
               D3.*Q, ...13
               D2.*T.^2, ...14
               D3.*T.^2, ...15
               D2.*Q.^2, ...16
               D3.*Q.^2, ...17
               D2.*T.*Q, ...18
               D3.*T.*Q, ...19
               (prestdz).*T, ...20    %%%%We use standardized pre-test score for numerical stability
               (prestdz).*T.^2, ...21 %%%%We use standardized pre-test score for numerical stability
               (prestdz).*Q, ...22    %%%%We use standardized pre-test score for numerical stability
               (prestdz).*Q.^2, ...23 %%%%We use standardized pre-test score for numerical stability
               (prestdz).*T.*Q] ; %24
labels_SR15  = {'Const*';
               'T';'T^2'; ...
               'Q';'Q^2';'T*Q'; ...
               'PreScore*';'D2';'D3'; 'D2*T'; 'D3*T'; 'D2*Q'; 'D3*Q'; 'D2*T^2'; 'D3*T^2'; 'D2*Q^2'; 'D3*Q^2'; 'D2*T*Q'; 'D3*T*Q'; ...
               'pre*T';'pre*T^2';'pre*Q';'pre*Q^2';'pre*T*Q'} ;  
[alpha_SR15,se_SR15,tstat_SR15,pval_SR15,CI_SR15,out_SR15] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,[],[],0,[],[],[]); 

StDevEffectATE    = nan(length(alpha_SR15)+2,1) ;
Delta0            = alpha_SR15(1) + prestdz*alpha_SR15(7) + D2*alpha_SR15(8) + D3*alpha_SR15(9) ; 
StDevEffectATE(1) = std(Delta0)/StDevDeltaScore  ;  %%%%This is a mean SD effect of Delta0

dDeltadT = (alpha_SR15(2) + D2*alpha_SR15(10) + D3*alpha_SR15(11) + prestdz*alpha_SR15(20)) ...
            + 2*(alpha_SR15(3) + D2*alpha_SR15(14) + D3*alpha_SR15(15) + prestdz*alpha_SR15(21)).*T ...
            + Q.*(alpha_SR15(6) + D2*alpha_SR15(18) + D3*alpha_SR15(19) + prestdz*alpha_SR15(24)) ; 
StDevEffectATE(2) = mean( BenchmarkChangeTstd*dDeltadT/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD partial effect of T

dDeltadQ = (alpha_SR15(4) + D2*alpha_SR15(12) + D3*alpha_SR15(13) + prestdz*alpha_SR15(22)) ...
            + 2*(alpha_SR15(5) + D2*alpha_SR15(16) + D3*alpha_SR15(17) + prestdz*alpha_SR15(23)).*Q ...
            + T.*(alpha_SR15(6) + D2*alpha_SR15(18) + D3*alpha_SR15(19) + prestdz*alpha_SR15(24)) ;
StDevEffectATE(4) = mean( BenchmarkChangeQstd*dDeltadQ/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD partial effect of Q   



StDevEffectATE(7) = mean( ( alpha_SR15(7) + T*alpha_SR15(20) + (T.^2)*alpha_SR15(21) + Q*alpha_SR15(22) + (Q.^2)*alpha_SR15(23) + T.*Q*alpha_SR15(24) - Presig )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
StDevEffectATE(8) = mean( ( alpha_SR15(8) + T*alpha_SR15(10) + Q*alpha_SR15(12) + (T.^2)*alpha_SR15(14) + (Q.^2)*alpha_SR15(16) + T.*Q*alpha_SR15(18) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(9) = mean( ( alpha_SR15(9) + T*alpha_SR15(11) + Q*alpha_SR15(13) + (T.^2)*alpha_SR15(15) + (Q.^2)*alpha_SR15(17) + T.*Q*alpha_SR15(19) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3


OtherEffect    = nan(length(alpha_SR15)+2,1) ;
%%%%This is a mean total derivative of kappa with respect to Q, for a
%%%%change of BenchmarkChangeQstd, as a fraction of StDevDeltaScore
OtherEffect(4) = mean( ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd/StDevDeltaScore ) ;  %%%%This is a mean total derivative effect of Q (in SD(DeltaS) units)
OtherEffect(7) = mean( ( T*alpha_SR15(20) + (T.^2)*alpha_SR15(21) + Q*alpha_SR15(22) + (Q.^2)*alpha_SR15(23) + T.*Q*alpha_SR15(24) )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score SLOPE EFFECTS ONLY (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
OtherEffect(8) = mean( ( T*alpha_SR15(10) + Q*alpha_SR15(12) + (T.^2)*alpha_SR15(14) + (Q.^2)*alpha_SR15(16) + T.*Q*alpha_SR15(18) )/StDevDeltaScore ) ;  %%%%This is an ATE change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(9) = mean( ( T*alpha_SR15(11) + Q*alpha_SR15(13) + (T.^2)*alpha_SR15(15) + (Q.^2)*alpha_SR15(17) + T.*Q*alpha_SR15(19) )/StDevDeltaScore ) ;  %%%%This is an ATE change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(10) = std( ( T*alpha_SR15(10) + Q*alpha_SR15(12) + (T.^2)*alpha_SR15(14) + (Q.^2)*alpha_SR15(16) + T.*Q*alpha_SR15(18) ) ) ;  %%%%This is a TE variance of a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(11) = std( ( T*alpha_SR15(11) + Q*alpha_SR15(13) + (T.^2)*alpha_SR15(15) + (Q.^2)*alpha_SR15(17) + T.*Q*alpha_SR15(19) ) ) ;  %%%%This is a TE variance of a change from District 1 to District 3 SLOPE EFFECTS ONLY
disp('-----------------------------------------------------------------------------------')
disp('-------------------------  SHORT-RUN PROGRESS: MODEL 1.5  -------------------------')
disp('----------  (THIS IS THE SPECIFICATION DISPLAYED IN COLUMN (1) OF TABLE 6)  -------')
disp('-----------------------------------------------------------------------------------')
baselinetemp = sort( alpha_SR15(1) + prestdz*alpha_SR15(7) + ...
                            + D2*alpha_SR15(8) + D3*alpha_SR15(9) - pre) ;
BASELINE_SR15 = [mean(baselinetemp) median(baselinetemp) std(baselinetemp) baselinetemp(round(length(baselinetemp)*0.1)) baselinetemp(round(length(baselinetemp)*0.9))] %#ok
out_SR15.results.StDevEffectATE = StDevEffectATE;
out_SR15.results.OtherEffect = OtherEffect;
WhiteTest  = out_SR15.White %#ok
Method     = out_SR15.method %#ok
disp(out_SR15.results(:,:)) 
%%%%Test for joint significance of all intercept terms:           
R = zeros(4,length(ExpVars_SR15(1,:))) ; R(1,1)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; r = zeros(4,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
InterceptTest = temp.Wald %#ok
%%%%Test for joint significance of all NON-S^{pre} intercept terms:           
R = zeros(3,length(ExpVars_SR15(1,:))) ; R(1,1)=1; R(2,8)=1; R(3,9)=1; r = zeros(3,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
InterceptTest_NonSpre = temp.Wald %#ok
%%%%Test for joint significance of all T terms:           
R = zeros(12,length(ExpVars_SR15(1,:))) ; R(1,2)=1; R(2,3)=1; R(3,6)=1; R(4,10)=1; R(5,11)=1; R(6,14)=1; R(7,15)=1; R(8,18)=1; R(9,19)=1; 
                                         R(10,20)=1; R(11,21)=1; R(12,24)=1;  r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
TTest = temp.Wald %#ok
%%%%Test for joint significance of all Q terms:           
R = zeros(12,length(ExpVars_SR15(1,:))) ; R(1,4)=1; R(2,5)=1; R(3,6)=1; R(4,12)=1; R(5,13)=1; R(6,16)=1; R(7,17)=1; R(8,18)=1; R(9,19)=1;
                                         R(10,22)=1; R(11,23)=1; R(12,24)=1;  r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
QTest = temp.Wald %#ok   
%%%%Test for joint significance of all District 2 terms:           
R = zeros(6,length(ExpVars_SR15(1,:))) ; R(1,8)=1; R(2,10)=1; R(3,12)=1; R(4,14)=1; R(5,16)=1; R(6,18)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
Dist2Test = temp.Wald %#ok
%%%%Test for joint significance of all District 2 slope interactions:           
R = zeros(5,length(ExpVars_SR15(1,:))) ; R(1,10)=1; R(2,12)=1; R(3,14)=1; R(4,16)=1; R(5,18)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
Dist2InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all District 3 terms:           
R = zeros(6,length(ExpVars_SR15(1,:))) ; R(1,9)=1; R(2,11)=1; R(3,13)=1; R(4,15)=1; R(5,17)=1; R(6,19)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
Dist3Test = temp.Wald %#ok
%%%%Test for joint significance of all District 3 slope interactions:           
R = zeros(5,length(ExpVars_SR15(1,:))) ; R(1,11)=1; R(2,13)=1; R(3,15)=1; R(4,17)=1; R(5,19)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
Dist3InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District terms:           
R = zeros(12,length(ExpVars_SR15(1,:))) ; for ii=1:12; R(ii,7+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
DistTest = temp.Wald %#ok
%%%%Test for joint significance of all School District slope interactions:           
R = zeros(10,length(ExpVars_SR15(1,:))) ; for ii=1:10; R(ii,9+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
DistInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District Intercept Terms:           
R = zeros(2,length(ExpVars_SR15(1,:))) ; R(1,8) = 1 ; R(2,9) = 1 ; r = zeros(2,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
DistInterceptTest = temp.Wald %#ok          
%%%%Test for joint significance of all pre-test terms:           
R = zeros(6,length(ExpVars_SR15(1,:))) ; R(1,7)=1; for ii=1:5; R(ii+1,19+ii)=1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
PreScoreTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test slope interactions:           
R = zeros(5,length(ExpVars_SR15(1,:))) ; for ii=1:5; R(ii,19+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR15,'FGLSrobust',labels_SR15,R,r,0,[],[],[]); 
PreScoreInteractionsTest = temp.Wald %#ok



%%


ExpVars_SR2 = [ones(size(Postscore)), ... 1
               T, ... 2
               T.^2, ... 3
               Q, ... 4
               Q.^2, ...5
               T.*Q, ...6
               prestdz, ...7 %%%%We use standardized pre-test score for numerical stability
               log(Y_EimputeShrunk), ...8
               log(Y_LimputeShrunk), ...9
               D2, ...10
               D3, ...11
               D2.*T, ...12
               D3.*T, ...13
               D2.*Q, ...14
               D3.*Q, ...15
               D2.*T.^2, ...16
               D3.*T.^2, ...17
               D2.*Q.^2, ...18
               D3.*Q.^2, ...19
               D2.*T.*Q, ...20
               D3.*T.*Q, ...21
               log(Y_EimputeShrunk).*T, ...22
               log(Y_LimputeShrunk).*T, ...23
               log(Y_EimputeShrunk).*Q, ...24
               log(Y_LimputeShrunk).*Q, ...25
               log(Y_EimputeShrunk).*T.^2, ...26
               log(Y_LimputeShrunk).*T.^2, ...27
               log(Y_EimputeShrunk).*Q.^2, ...28
               log(Y_LimputeShrunk).*Q.^2, ...29
               log(Y_EimputeShrunk).*T.*Q, ...30
               log(Y_LimputeShrunk).*T.*Q, ...31
               (prestdz).*T, ...32    %%%%We use standardized pre-test score for numerical stability
               (prestdz).*T.^2, ...33 %%%%We use standardized pre-test score for numerical stability
               (prestdz).*Q, ...34    %%%%We use standardized pre-test score for numerical stability
               (prestdz).*Q.^2, ...35 %%%%We use standardized pre-test score for numerical stability
               (prestdz).*T.*Q] ; %36
labels_SR2  = {'Const*';
               'T';'T^2'; ...
               'Q';'Q^2';'T*Q'; ...
               'PreScore*';'log(thetaE)';'log(thetaL)'; ...
               'D2';'D3'; 'D2*T'; 'D3*T'; 'D2*Q'; 'D3*Q'; 'D2*T^2'; 'D3*T^2'; 'D2*Q^2'; 'D3*Q^2'; 'D2*T*Q'; 'D3*T*Q'; ...
               'log(thetaE)*T'; 'log(thetaL)*T'; 'log(thetaE)*Q'; 'log(thetaL)*Q'; 'log(thetaE)*T^2'; 'log(thetaL)*T^2'; 'log(thetaE)*Q^2'; 'log(thetaL)*Q^2'; 'log(thetaE)*T*Q'; 'log(thetaL)*T*Q'; ...
               'pre*T';'pre*T^2';'pre*Q';'pre*Q^2';'pre*T*Q'} ;
[alpha_SR2,se_SR2,tstat_SR2,pval_SR2,CI_SR2,out_SR2] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,[],[],0,[],[],[]); 

StDevEffectATE    = nan(length(alpha_SR2)+2,1) ;
Delta0            = alpha_SR2(1) + log(Y_EimputeShrunk)*alpha_SR2(8) + ...
                    log(Y_LimputeShrunk)*alpha_SR2(9) + D2*alpha_SR2(10) + D3*alpha_SR2(11) ;
StDevEffectATE(1) = std(Delta0)/StDevDeltaScore  ;  %%%%This is a mean SD effect of Delta0

dDeltadT = (alpha_SR2(2) + D2*alpha_SR2(12) + D3*alpha_SR2(13) + log(Y_EimputeShrunk)*alpha_SR2(22) + log(Y_LimputeShrunk)*alpha_SR2(23) + prestdz*alpha_SR2(32)) ...
            + 2*(alpha_SR2(3) + D2*alpha_SR2(16) + D3*alpha_SR2(17) + log(Y_EimputeShrunk)*alpha_SR2(26) + log(Y_LimputeShrunk)*alpha_SR2(27) + prestdz*alpha_SR2(33)).*T ...
            + Q.*(alpha_SR2(6) + D2*alpha_SR2(20) + D3*alpha_SR2(21) + log(Y_EimputeShrunk)*alpha_SR2(30) + log(Y_LimputeShrunk)*alpha_SR2(31) + prestdz*alpha_SR2(36)) ;
StDevEffectATE(2) = mean( BenchmarkChangeTstd*dDeltadT/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD partial effect of T

dDeltadQ = (alpha_SR2(4) + D2*alpha_SR2(14) + D3*alpha_SR2(15) + log(Y_EimputeShrunk)*alpha_SR2(24) + log(Y_LimputeShrunk)*alpha_SR2(25) + prestdz*alpha_SR2(34)) ...
            + 2*(alpha_SR2(5) + D2*alpha_SR2(18) + D3*alpha_SR2(19) + log(Y_EimputeShrunk)*alpha_SR2(28) + log(Y_LimputeShrunk)*alpha_SR2(29) + prestdz*alpha_SR2(35)).*Q ...
            + T.*(alpha_SR2(6) + D2*alpha_SR2(20) + D3*alpha_SR2(21) + log(Y_EimputeShrunk)*alpha_SR2(30) + log(Y_LimputeShrunk)*alpha_SR2(31) + prestdz*alpha_SR2(36)) ;
StDevEffectATE(4) = mean( BenchmarkChangeQstd*dDeltadQ/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD partial effect of Q   



StDevEffectATE(7) = mean( ( alpha_SR2(7) + T*alpha_SR2(32) + (T.^2)*alpha_SR2(33) + Q*alpha_SR2(34) + (Q.^2)*alpha_SR2(35) + T.*Q*alpha_SR2(36) - Presig )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
StDevEffectATE(8) = mean( std(log(Y_EimputeShrunk))*( alpha_SR2(8) + T*alpha_SR2(22) + Q*alpha_SR2(24) + (T.^2)*alpha_SR2(26)+ (Q.^2)*alpha_SR2(28) + T.*Q*alpha_SR2(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE)
StDevEffectATE(9) = mean( std(log(Y_LimputeShrunk))*( alpha_SR2(9) + T*alpha_SR2(23) + Q*alpha_SR2(25) + (T.^2)*alpha_SR2(27)+ (Q.^2)*alpha_SR2(29) + T.*Q*alpha_SR2(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL)
StDevEffectATE(10) = mean( ( alpha_SR2(10) + T*alpha_SR2(12) + Q*alpha_SR2(14) + (T.^2)*alpha_SR2(16) + (Q.^2)*alpha_SR2(18) + T.*Q*alpha_SR2(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(11) = mean( ( alpha_SR2(11) + T*alpha_SR2(13) + Q*alpha_SR2(15) + (T.^2)*alpha_SR2(17) + (Q.^2)*alpha_SR2(19) + T.*Q*alpha_SR2(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3


OtherEffect    = nan(length(alpha_SR2)+2,1) ;
%%%%This is a mean total derivative of kappa with respect to Q, for a
%%%%change of BenchmarkChangeQstd, as a fraction of StDevDeltaScore
OtherEffect(4) = mean( ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd/StDevDeltaScore ) ;  %%%%This is a mean total derivative effect of Q (in original DeltaS units)
OtherEffect(7) = mean( ( T*alpha_SR2(32) + (T.^2)*alpha_SR2(33) + Q*alpha_SR2(34) + (Q.^2)*alpha_SR2(35) + T.*Q*alpha_SR2(36) )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score SLOPE EFFECTS ONLY (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
OtherEffect(8) = mean( std(log(Y_EimputeShrunk))*( T*alpha_SR2(22) + Q*alpha_SR2(24) + (T.^2)*alpha_SR2(26)+ (Q.^2)*alpha_SR2(28) + T.*Q*alpha_SR2(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE) SLOPE EFFECTS ONLY
OtherEffect(9) = mean( std(log(Y_LimputeShrunk))*( T*alpha_SR2(23) + Q*alpha_SR2(25) + (T.^2)*alpha_SR2(27)+ (Q.^2)*alpha_SR2(29) + T.*Q*alpha_SR2(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL) SLOPE EFFECTS ONLY
OtherEffect(10) = mean( ( T*alpha_SR2(12) + Q*alpha_SR2(14) + (T.^2)*alpha_SR2(16) + (Q.^2)*alpha_SR2(18) + T.*Q*alpha_SR2(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(11) = mean( ( T*alpha_SR2(13) + Q*alpha_SR2(15) + (T.^2)*alpha_SR2(17) + (Q.^2)*alpha_SR2(19) + T.*Q*alpha_SR2(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(12) = std( ( T*alpha_SR2(12) + Q*alpha_SR2(14) + (T.^2)*alpha_SR2(16) + (Q.^2)*alpha_SR2(18) + T.*Q*alpha_SR2(20) ) ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(13) = std( ( T*alpha_SR2(13) + Q*alpha_SR2(15) + (T.^2)*alpha_SR2(17) + (Q.^2)*alpha_SR2(19) + T.*Q*alpha_SR2(21) ) ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
disp('-----------------------------------------------------------------------------------')
disp('-------------------------  SHORT-RUN PROGRESS: MODEL 2  ---------------------------')
disp('-----------------------------------------------------------------------------------')
baselinetemp = sort( alpha_SR2(1)+prestdz*alpha_SR2(7)+log(Y_EimputeShrunk)*alpha_SR2(8)+log(Y_LimputeShrunk)*alpha_SR2(9) ...
                            +D2*alpha_SR2(10)+D3*alpha_SR2(11) - pre) ;
BASELINE_SR2 = [mean(baselinetemp) median(baselinetemp) std(baselinetemp) baselinetemp(round(length(baselinetemp)*0.1)) baselinetemp(round(length(baselinetemp)*0.9))] %#ok
out_SR2.results.StDevEffectATE = StDevEffectATE;
out_SR2.results.OtherEffect = OtherEffect;
WhiteTest  = out_SR2.White %#ok
Method     = out_SR2.method %#ok
disp(out_SR2.results(:,:)) 
%%%%Test for joint significance of all intercept terms:           
R = zeros(6,length(ExpVars_SR2(1,:))) ; R(1,1)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; R(5,10)=1; R(6,11)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
InterceptTest = temp.Wald %#ok
%%%%Test for joint significance of all NON-S^{pre} intercept terms:           
R = zeros(5,length(ExpVars_SR2(1,:))) ; R(1,1)=1; R(2,8)=1; R(3,9)=1; R(4,10)=1; R(5,11)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
InterceptTest_NonSpre = temp.Wald %#ok
%%%%Test for joint significance of all T terms:           
R = zeros(18,length(ExpVars_SR2(1,:))) ; R(1,2)=1; R(2,3)=1; R(3,6)=1; R(4,12)=1; R(5,13)=1; R(6,16)=1; R(7,17)=1; R(8,20)=1; R(9,21)=1; 
                                         R(10,22)=1; R(11,23)=1; R(12,26)=1; R(13,27)=1; R(14,30)=1; R(15,31)=1; R(16,32)=1; R(17,33)=1; R(18,36)=1;  r = zeros(18,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
TTest = temp.Wald %#ok
%%%%Test for joint significance of all Q terms:           
R = zeros(18,length(ExpVars_SR2(1,:))) ; R(1,4)=1; R(2,5)=1; R(3,6)=1; R(4,14)=1; R(5,15)=1; R(6,18)=1; R(7,19)=1; R(8,20)=1; R(9,21)=1;
                                         R(10,24)=1; R(11,25)=1; R(12,28)=1; R(13,29)=1; R(14,30)=1; R(15,31)=1; R(16,34)=1; R(17,35)=1; R(18,36)=1;  r = zeros(18,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
QTest = temp.Wald %#ok   
%%%%Test for joint significance of all District 2 terms:           
R = zeros(6,length(ExpVars_SR2(1,:))) ; R(1,10)=1; R(2,12)=1; R(3,14)=1; R(4,16)=1; R(5,18)=1; R(6,20)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
Dist2Test = temp.Wald %#ok
%%%%Test for joint significance of all District 2 slope interactions:           
R = zeros(5,length(ExpVars_SR2(1,:))) ; R(1,12)=1; R(2,14)=1; R(3,16)=1; R(4,18)=1; R(5,20)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
Dist2InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all District 3 terms:           
R = zeros(6,length(ExpVars_SR2(1,:))) ; R(1,11)=1; R(2,13)=1; R(3,15)=1; R(4,17)=1; R(5,19)=1; R(6,21)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
Dist3Test = temp.Wald %#ok
%%%%Test for joint significance of all District 3 slope interactions:           
R = zeros(5,length(ExpVars_SR2(1,:))) ; R(1,13)=1; R(2,15)=1; R(3,17)=1; R(4,19)=1; R(5,21)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
Dist3InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District terms:           
R = zeros(12,length(ExpVars_SR2(1,:))) ; for ii=1:12; R(ii,9+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
DistTest = temp.Wald %#ok
%%%%Test for joint significance of all School District slope interactions:           
R = zeros(10,length(ExpVars_SR2(1,:))) ; for ii=1:10; R(ii,11+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
DistInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District Intercept Terms:           
R = zeros(2,length(ExpVars_SR2(1,:))) ; R(1,10) = 1 ; R(2,11) = 1 ; r = zeros(2,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
DistInterceptTest = temp.Wald %#ok          
%%%%Test for joint significance of all ThetaE terms: 
R = zeros(6,length(ExpVars_SR2(1,:))) ; R(1,8)=1 ; R(2,22)=1 ; R(3,24)=1 ; R(4,26)=1 ; R(5,28)=1 ; R(6,30)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
ThetaETest = temp.Wald %#ok           
%%%%Test for joint significance of all ThetaE interaction terms: 
R = zeros(5,length(ExpVars_SR2(1,:))) ; R(1,22)=1 ; R(2,24)=1 ; R(3,26)=1 ; R(4,28)=1 ; R(5,30)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
ThetaEInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL terms:          
R = zeros(6,length(ExpVars_SR2(1,:))) ; R(1,9)=1 ; R(2,23)=1 ; R(3,25)=1 ; R(4,27)=1 ; R(5,29)=1 ; R(6,31)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
ThetaLTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL interaction terms:          
R = zeros(5,length(ExpVars_SR2(1,:))) ; R(1,23)=1 ; R(2,25)=1 ; R(3,27)=1 ; R(4,29)=1 ; R(5,31)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
ThetaLInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all Student Type slope interactions:          
R = zeros(10,length(ExpVars_SR2(1,:))) ; for ii=1:10; R(ii,21+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
TypeInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test terms:           
R = zeros(6,length(ExpVars_SR2(1,:))) ; R(1,7)=1; for ii=1:5; R(ii+1,31+ii)=1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
PreScoreTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test slope interactions:           
R = zeros(5,length(ExpVars_SR2(1,:))) ; for ii=1:5; R(ii,31+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR2,'FGLSrobust',labels_SR2,R,r,0,[],[],[]); 
PreScoreInteractionsTest = temp.Wald %#ok

%%




ExpVars_SR3 = [ones(size(Postscore)), ... 1
               T, ... 2
               T.^2, ... 3
               Q, ... 4
               Q.^2, ...5
               T.*Q, ...6
               prestdz, ...7 %%%%We use standardized pre-test score for numerical stability
               log(Y_EimputeShrunk), ...8
               log(Y_LimputeShrunk), ...9
               D2, ...10
               D3, ...11
               D2.*T, ...12
               D3.*T, ...13
               D2.*Q, ...14
               D3.*Q, ...15
               D2.*T.^2, ...16
               D3.*T.^2, ...17
               D2.*Q.^2, ...18
               D3.*Q.^2, ...19
               D2.*T.*Q, ...20
               D3.*T.*Q, ...21
               log(Y_EimputeShrunk).*T, ...22
               log(Y_LimputeShrunk).*T, ...23
               log(Y_EimputeShrunk).*Q, ...24
               log(Y_LimputeShrunk).*Q, ...25
               log(Y_EimputeShrunk).*T.^2, ...26
               log(Y_LimputeShrunk).*T.^2, ...27
               log(Y_EimputeShrunk).*Q.^2, ...28
               log(Y_LimputeShrunk).*Q.^2, ...29
               log(Y_EimputeShrunk).*T.*Q, ...30
               log(Y_LimputeShrunk).*T.*Q, ...31
               (prestdz).*T, ...32    %%%%We use standardized pre-test score for numerical stability
               (prestdz).*T.^2, ...33 %%%%We use standardized pre-test score for numerical stability
               (prestdz).*Q, ...34    %%%%We use standardized pre-test score for numerical stability
               (prestdz).*Q.^2, ...35 %%%%We use standardized pre-test score for numerical stability
               (prestdz).*T.*Q, ...36 %%%%We use standardized pre-test score for numerical stability
               grade5, ...37
               grade5.*T, ...38
               grade5.*T.^2, ...39
               grade5.*Q, ...40
               grade5.*Q.^2, ...41
               grade5.*T.*Q, ... 42
               fem, ...43
               fem.*T, ...44
               fem.*T.^2, ...45
               fem.*Q, ...46
               fem.*Q.^2, ...47
               fem.*T.*Q, ...48
               blk, ...49
               blk.*T, ...50
               blk.*T.^2, ...51
               blk.*Q, ...52
               blk.*Q.^2, ...53
               blk.*T.*Q, ...54
               hsp, ...55
               hsp.*T, ...56
               hsp.*T.^2, ...57
               hsp.*Q, ...58
               hsp.*Q.^2, ...59
               hsp.*T.*Q] ; %60
labels_SR3  = {'Const*';
               'T';'T^2'; ...
               'Q';'Q^2';'T*Q'; ...
               'PreScore*';'log(thetaE)';'log(thetaL)'; ...
               'D2';'D3'; 'D2*T'; 'D3*T'; 'D2*Q'; 'D3*Q'; 'D2*T^2'; 'D3*T^2'; 'D2*Q^2'; 'D3*Q^2'; 'D2*T*Q'; 'D3*T*Q'; ...
               'log(thetaE)*T'; 'log(thetaL)*T'; 'log(thetaE)*Q'; 'log(thetaL)*Q'; 'log(thetaE)*T^2'; 'log(thetaL)*T^2'; 'log(thetaE)*Q^2'; 'log(thetaL)*Q^2'; 'log(thetaE)*T*Q'; 'log(thetaL)*T*Q'; ...
               'pre*T';'pre*T^2';'pre*Q';'pre*Q^2';'pre*T*Q';'grade5';'grade5*T';'grade5*T^2';'grade5*Q';'grade5*Q^2';'grade5*T*Q';...
               'fem';'fem*T';'fem*T^2';'fem*Q';'fem*Q^2';'fem*T*Q';'blk';'blk*T';'blk*T^2';'blk*Q';'blk*Q^2';'blk*T*Q';'hsp';'hsp*T';'hsp*T^2';'hsp*Q';'hsp*Q^2';'hsp*T*Q'} ;
[alpha_SR3,se_SR3,tstat_SR3,pval_SR3,CI_SR3,out_SR3] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,[],[],0,[],[],[]); 

StDevEffectATE    = nan(length(alpha_SR3)+2,1) ;
Delta0 = alpha_SR3(1) + log(Y_EimputeShrunk)*alpha_SR3(8) + log(Y_LimputeShrunk)*alpha_SR3(9) ...
         + D2*alpha_SR3(10) + D3*alpha_SR3(11) + grade5*alpha_SR3(37) + fem*alpha_SR3(43) + blk*alpha_SR3(49) + hsp*alpha_SR3(55) ;
StDevEffectATE(1) = std(Delta0)/StDevDeltaScore  ;

dDeltadT = (alpha_SR3(2) + D2*alpha_SR3(12) + D3*alpha_SR3(13) + log(Y_EimputeShrunk)*alpha_SR3(22) + log(Y_LimputeShrunk)*alpha_SR3(23) ...
              + prestdz*alpha_SR3(32) + grade5*alpha_SR3(38) + fem*alpha_SR3(44) + blk*alpha_SR3(50) + hsp*alpha_SR3(56)) ...
            + 2*(alpha_SR3(3) + D2*alpha_SR3(16) + D3*alpha_SR3(17) + log(Y_EimputeShrunk)*alpha_SR3(26) + log(Y_LimputeShrunk)*alpha_SR3(27) ...
              + prestdz*alpha_SR3(33) + grade5*alpha_SR3(39) + fem*alpha_SR3(45) + blk*alpha_SR3(51) + hsp*alpha_SR3(57)).*T ...
            + Q.*(alpha_SR3(6) + D2*alpha_SR3(20) + D3*alpha_SR3(21) + log(Y_EimputeShrunk)*alpha_SR3(30) + log(Y_LimputeShrunk)*alpha_SR3(31) ...
              + prestdz*alpha_SR3(36) + grade5*alpha_SR3(42) + fem*alpha_SR3(48) + blk*alpha_SR3(54) + hsp*alpha_SR3(60)) ;
StDevEffectATE(2) = mean( BenchmarkChangeTstd*dDeltadT/StDevDeltaScore ) ;  %%%%This is mean pseudo-SD effect of T

dDeltadQ = (alpha_SR3(4) + D2*alpha_SR3(14) + D3*alpha_SR3(15) + log(Y_EimputeShrunk)*alpha_SR3(24) + log(Y_LimputeShrunk)*alpha_SR3(25) ...
              + prestdz*alpha_SR3(34) + grade5*alpha_SR3(40) + fem*alpha_SR3(46) + blk*alpha_SR3(52) + hsp*alpha_SR3(58)) ...
            + 2*(alpha_SR3(5) + D2*alpha_SR3(18) + D3*alpha_SR3(19) + log(Y_EimputeShrunk)*alpha_SR3(28) + log(Y_LimputeShrunk)*alpha_SR3(29) ...
              + prestdz*alpha_SR3(35) + grade5*alpha_SR3(41) + fem*alpha_SR3(47) + blk*alpha_SR3(53) + hsp*alpha_SR3(59)).*Q ...
            + T.*(alpha_SR3(6) + D2*alpha_SR3(20) + D3*alpha_SR3(21) + log(Y_EimputeShrunk)*alpha_SR3(30) + log(Y_LimputeShrunk)*alpha_SR3(31) ...
              + prestdz*alpha_SR3(36) + grade5*alpha_SR3(42) + fem*alpha_SR3(48) + blk*alpha_SR3(54) + hsp*alpha_SR3(60)) ;
StDevEffectATE(4) = mean( BenchmarkChangeQstd*dDeltadQ/StDevDeltaScore ) ;  %%%%This is a standard deviation change in Q

StDevEffectATE(7) = mean( ( alpha_SR3(7) + T*alpha_SR3(32) + (T.^2)*alpha_SR3(33) + Q*alpha_SR3(34) + (Q.^2)*alpha_SR3(35) + T.*Q*alpha_SR3(36) - Presig )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
StDevEffectATE(8) = mean( std(log(Y_EimputeShrunk))*( alpha_SR3(8) + T*alpha_SR3(22) + Q*alpha_SR3(24) + (T.^2)*alpha_SR3(26)+ (Q.^2)*alpha_SR3(28) + T.*Q*alpha_SR3(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE)
StDevEffectATE(9) = mean( std(log(Y_LimputeShrunk))*( alpha_SR3(9) + T*alpha_SR3(23) + Q*alpha_SR3(25) + (T.^2)*alpha_SR3(27)+ (Q.^2)*alpha_SR3(29) + T.*Q*alpha_SR3(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL)
StDevEffectATE(10) = mean( ( alpha_SR3(10) + T*alpha_SR3(12) + Q*alpha_SR3(14) + (T.^2)*alpha_SR3(16) + (Q.^2)*alpha_SR3(18) + T.*Q*alpha_SR3(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(11) = mean( ( alpha_SR3(11) + T*alpha_SR3(13) + Q*alpha_SR3(15) + (T.^2)*alpha_SR3(17) + (Q.^2)*alpha_SR3(19) + T.*Q*alpha_SR3(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3
StDevEffectATE(37) = mean( ( alpha_SR3(37) + T*alpha_SR3(38) + (T.^2)*alpha_SR3(39) + Q*alpha_SR3(40) + (Q.^2)*alpha_SR3(41) + T.*Q*alpha_SR3(42) )/StDevDeltaScore ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectATE(43) = mean( ( alpha_SR3(43) + T*alpha_SR3(44) + (T.^2)*alpha_SR3(45) + Q*alpha_SR3(46) + (Q.^2)*alpha_SR3(47) + T.*Q*alpha_SR3(48) )/StDevDeltaScore ) ;  %%%%This is a change from male to female
StDevEffectATE(49) = mean( ( alpha_SR3(49) + T*alpha_SR3(50) + (T.^2)*alpha_SR3(51) + Q*alpha_SR3(52) + (Q.^2)*alpha_SR3(53) + T.*Q*alpha_SR3(54) )/StDevDeltaScore ) ;  %%%%This is a change from wao to black
StDevEffectATE(55) = mean( ( alpha_SR3(55) + T*alpha_SR3(56) + (T.^2)*alpha_SR3(57) + Q*alpha_SR3(58) + (Q.^2)*alpha_SR3(59) + T.*Q*alpha_SR3(60) )/StDevDeltaScore ) ;  %%%%This is a change from wao to hispanic

OtherEffect    = nan(length(alpha_SR3)+2,1) ;
%%%%This is a mean total derivative of kappa with respect to Q, for a
%%%%change of BenchmarkChangeQstd, as a fraction of StDevDeltaScore
OtherEffect(4) = mean( ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd/StDevDeltaScore ) ;  %%%%This is a mean total derivative effect of Q (in original DeltaS units)
OtherEffect(7) = mean( ( T*alpha_SR3(32) + (T.^2)*alpha_SR3(33) + Q*alpha_SR3(34) + (Q.^2)*alpha_SR3(35) + T.*Q*alpha_SR3(36) )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score SLOPE EFFECTS ONLY (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
OtherEffect(8) = mean( std(log(Y_EimputeShrunk))*( T*alpha_SR3(22) + Q*alpha_SR3(24) + (T.^2)*alpha_SR3(26)+ (Q.^2)*alpha_SR3(28) + T.*Q*alpha_SR3(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE) SLOPE EFFECTS ONLY
OtherEffect(9) = mean( std(log(Y_LimputeShrunk))*( T*alpha_SR3(23) + Q*alpha_SR3(25) + (T.^2)*alpha_SR3(27)+ (Q.^2)*alpha_SR3(29) + T.*Q*alpha_SR3(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL) SLOPE EFFECTS ONLY
OtherEffect(10) = mean( ( T*alpha_SR3(12) + Q*alpha_SR3(14) + (T.^2)*alpha_SR3(16) + (Q.^2)*alpha_SR3(18) + T.*Q*alpha_SR3(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(11) = mean( ( T*alpha_SR3(13) + Q*alpha_SR3(15) + (T.^2)*alpha_SR3(17) + (Q.^2)*alpha_SR3(19) + T.*Q*alpha_SR3(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(12) = std( ( T*alpha_SR3(12) + Q*alpha_SR3(14) + (T.^2)*alpha_SR3(16) + (Q.^2)*alpha_SR3(18) + T.*Q*alpha_SR3(20) ) ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(13) = std( ( T*alpha_SR3(13) + Q*alpha_SR3(15) + (T.^2)*alpha_SR3(17) + (Q.^2)*alpha_SR3(19) + T.*Q*alpha_SR3(21) ) ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(37) = mean( ( T*alpha_SR3(38) + (T.^2)*alpha_SR3(39) + Q*alpha_SR3(40) + (Q.^2)*alpha_SR3(41) + T.*Q*alpha_SR3(42) )/StDevDeltaScore ) ;  %%%%This is a change from grade 6 to grade 5 SLOPE EFFECTS ONLY
disp('-----------------------------------------------------------------------------------')
disp('-------------------------  SHORT-RUN PROGRESS: MODEL 3  ---------------------------')
disp('-----------------------------------------------------------------------------------')
baselinetemp = sort( alpha_SR3(1)+prestdz*alpha_SR3(7)+log(Y_EimputeShrunk)*alpha_SR3(8)+log(Y_LimputeShrunk)*alpha_SR3(9) ...
                            +D2*alpha_SR3(10)+D3*alpha_SR3(11)+grade5*alpha_SR3(37)+fem*alpha_SR3(43)+blk*alpha_SR3(49)+hsp*alpha_SR3(55) - pre ) ;
BASELINE_SR3 = [mean(baselinetemp) median(baselinetemp) std(baselinetemp) baselinetemp(round(length(baselinetemp)*0.1)) baselinetemp(round(length(baselinetemp)*0.9))] %#ok
out_SR3.results.StDevEffectATE = StDevEffectATE;
out_SR3.results.OtherEffect = OtherEffect;
WhiteTest  = out_SR3.White %#ok
Method     = out_SR3.method %#ok
disp(out_SR3.results(:,:)) 
%%%%Test for joint significance of all intercept terms:           
R = zeros(10,length(ExpVars_SR3(1,:))) ; R(1,1)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; R(5,10)=1; R(6,11)=1; R(7,37)=1; R(8,43)=1; R(9,49)=1; R(10,55)=1; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
InterceptTest = temp.Wald %#ok
%%%%Test for joint significance of all NON-S^{pre} intercept terms:           
R = zeros(9,length(ExpVars_SR3(1,:))) ; R(1,1)=1; R(2,8)=1; R(3,9)=1; R(4,10)=1; R(5,11)=1; R(6,37)=1; R(7,43)=1; R(8,49)=1; R(9,55)=1; r = zeros(9,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
InterceptTest_NonSpre = temp.Wald %#ok
%%%%Test for joint significance of all T terms:           
R = zeros(30,length(ExpVars_SR3(1,:))) ; R(1,2)=1; R(2,3)=1; R(3,6)=1; R(4,12)=1; R(5,13)=1; R(6,16)=1; R(7,17)=1; R(8,20)=1; R(9,21)=1; 
                                         R(10,22)=1; R(11,23)=1; R(12,26)=1; R(13,27)=1; R(14,30)=1; R(15,31)=1; R(16,32)=1; R(17,33)=1; R(18,36)=1;
                                         R(19,38)=1; R(20,39)=1; R(21,42)=1; R(22,44)=1; R(23,45)=1; R(24,48)=1;
                                         R(25,50)=1; R(26,51)=1; R(27,54)=1; R(28,56)=1; R(29,57)=1; R(30,60)=1;  r = zeros(30,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
TTest = temp.Wald %#ok
%%%%Test for joint significance of all Q terms:           
R = zeros(30,length(ExpVars_SR3(1,:))) ; R(1,4)=1; R(2,5)=1; R(3,6)=1; R(4,14)=1; R(5,15)=1; R(6,18)=1; R(7,19)=1; R(8,20)=1; R(9,21)=1;
                                         R(10,24)=1; R(11,25)=1; R(12,28)=1; R(13,29)=1; R(14,30)=1; R(15,31)=1; R(16,34)=1; R(17,35)=1; R(18,36)=1;
                                         R(19,40)=1; R(20,41)=1; R(21,42)=1; R(22,46)=1; R(23,47)=1; R(24,48)=1;
                                         R(25,52)=1; R(26,53)=1; R(27,54)=1; R(28,58)=1; R(29,59)=1; R(30,60)=1;  r = zeros(30,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
QTest = temp.Wald %#ok   
%%%%Test for joint significance of all District 2 terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; R(1,10)=1; R(2,12)=1; R(3,14)=1; R(4,16)=1; R(5,18)=1; R(6,20)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
Dist2Test = temp.Wald %#ok
%%%%Test for joint significance of all District 2 slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; R(1,12)=1; R(2,14)=1; R(3,16)=1; R(4,18)=1; R(5,20)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
Dist2InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all District 3 terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; R(1,11)=1; R(2,13)=1; R(3,15)=1; R(4,17)=1; R(5,19)=1; R(6,21)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
Dist3Test = temp.Wald %#ok
%%%%Test for joint significance of all District 3 slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; R(1,13)=1; R(2,15)=1; R(3,17)=1; R(4,19)=1; R(5,21)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
Dist3InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District terms:           
R = zeros(12,length(ExpVars_SR3(1,:))) ; for ii=1:12; R(ii,9+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
DistTest = temp.Wald %#ok
%%%%Test for joint significance of all School District slope interactions:           
R = zeros(10,length(ExpVars_SR3(1,:))) ; for ii=1:10; R(ii,11+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
DistInteractionsTest = temp.Wald %#ok    
%%%%Test for joint significance of all School District Intercept Terms:           
R = zeros(2,length(ExpVars_SR3(1,:))) ; R(1,10) = 1 ; R(2,11) = 1 ; r = zeros(2,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
DistInterceptTest = temp.Wald %#ok        
%%%%Test for joint significance of all ThetaE terms: 
R = zeros(6,length(ExpVars_SR3(1,:))) ; R(1,8)=1 ; R(2,22)=1 ; R(3,24)=1 ; R(4,26)=1 ; R(5,28)=1 ; R(6,30)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
ThetaETest = temp.Wald %#ok           
%%%%Test for joint significance of all ThetaE interaction terms: 
R = zeros(5,length(ExpVars_SR3(1,:))) ; R(1,22)=1 ; R(2,24)=1 ; R(3,26)=1 ; R(4,28)=1 ; R(5,30)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
ThetaEInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL terms:          
R = zeros(6,length(ExpVars_SR3(1,:))) ; R(1,9)=1 ; R(2,23)=1 ; R(3,25)=1 ; R(4,27)=1 ; R(5,29)=1 ; R(6,31)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
ThetaLTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL interaction terms:          
R = zeros(5,length(ExpVars_SR3(1,:))) ; R(1,23)=1 ; R(2,25)=1 ; R(3,27)=1 ; R(4,29)=1 ; R(5,31)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
ThetaLInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all Student Type slope interactions:          
R = zeros(10,length(ExpVars_SR3(1,:))) ; for ii=1:10; R(ii,21+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
TypeInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; R(1,7)=1; for ii=1:5; R(ii+1,31+ii)=1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
PreScoreTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; for ii=1:5; R(ii,31+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
PreScoreInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all age-cohort terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; for ii=1:6; R(ii,36+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
Grade5Test = temp.Wald %#ok
%%%%Test for joint significance of all age-cohort slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; for ii=1:5; R(ii,37+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
Grade5InteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all fem terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; for ii=1:6; R(ii,42+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
FemTest = temp.Wald %#ok
%%%%Test for joint significance of all fem slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; for ii=1:5; R(ii,43+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
FemInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all blk terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; for ii=1:6; R(ii,48+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
BlkTest = temp.Wald %#ok
%%%%Test for joint significance of all blk slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; for ii=1:5; R(ii,49+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
BlkInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all hsp terms:           
R = zeros(6,length(ExpVars_SR3(1,:))) ; for ii=1:6; R(ii,54+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
HspTest = temp.Wald %#ok
%%%%Test for joint significance of all hsp slope interactions:           
R = zeros(5,length(ExpVars_SR3(1,:))) ; for ii=1:5; R(ii,55+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR3,'FGLSrobust',labels_SR3,R,r,0,[],[],[]); 
HspInteractionsTest = temp.Wald %#ok
%%





ExpVars_SR4 = [ones(size(Postscore)), ... 1
               T, ... 2
               T.^2, ... 3
               Q, ... 4
               Q.^2, ...5
               T.*Q, ...6
               prestdz, ...7 %%%%We use standardized pre-test score for numerical stability in Model 6
               log(Y_EimputeShrunk), ...8
               log(Y_LimputeShrunk), ...9
               D2, ...10
               D3, ...11
               D2.*T, ...12
               D3.*T, ...13
               D2.*Q, ...14
               D3.*Q, ...15
               D2.*T.^2, ...16
               D3.*T.^2, ...17
               D2.*Q.^2, ...18
               D3.*Q.^2, ...19
               D2.*T.*Q, ...20
               D3.*T.*Q, ...21
               log(Y_EimputeShrunk).*T, ...22
               log(Y_LimputeShrunk).*T, ...23
               log(Y_EimputeShrunk).*Q, ...24
               log(Y_LimputeShrunk).*Q, ...25
               log(Y_EimputeShrunk).*T.^2, ...26
               log(Y_LimputeShrunk).*T.^2, ...27
               log(Y_EimputeShrunk).*Q.^2, ...28
               log(Y_LimputeShrunk).*Q.^2, ...29
               log(Y_EimputeShrunk).*T.*Q, ...30
               log(Y_LimputeShrunk).*T.*Q, ...31
               (prestdz).*T, ...32    %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*T.^2, ...33 %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*Q, ...34    %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*Q.^2, ...35 %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*T.*Q, ...36 %%%%We use standardized pre-test score for numerical stability in Model 6
               grade5, ...37
               grade5.*T, ...38
               grade5.*T.^2, ...39
               grade5.*Q, ...40
               grade5.*Q.^2, ...41
               grade5.*T.*Q, ... 42
               fem, ...43
               fem.*T, ...44
               fem.*T.^2, ...45
               fem.*Q, ...46
               fem.*Q.^2, ...47
               fem.*T.*Q, ...48
               blk, ...49
               blk.*T, ...50
               blk.*T.^2, ...51
               blk.*Q, ...52
               blk.*Q.^2, ...53
               blk.*T.*Q, ...54
               hsp, ...55
               hsp.*T, ...56
               hsp.*T.^2, ...57
               hsp.*Q, ...58
               hsp.*Q.^2, ...59
               hsp.*T.*Q, ...60
               SSEInc, ...61
               SSEInc.*T, ...62
               SSEInc.*T.^2, ...63
               SSEInc.*Q, ...64
               SSEInc.*Q.^2, ...65
               SSEInc.*T.*Q, ...66
               SSEIns, ...67
               SSEIns.*T, ...68
               SSEIns.*T.^2, ...69
               SSEIns.*Q, ...70
               SSEIns.*Q.^2, ...71
               SSEIns.*T.*Q, ...72
               NHelpers, ...73
               NHelpers.*T, ...74
               NHelpers.*T.^2, ...75
               NHelpers.*Q, ...76
               NHelpers.*Q.^2, ...77
               NHelpers.*T.*Q, ... 78
               NoHmIntrnt, ... 79
               NoHmIntrnt.*T, ... 80
               NoHmIntrnt.*T.^2, ... 81
               NoHmIntrnt.*Q, ... 82
               NoHmIntrnt.*Q.^2, ... 83
               NoHmIntrnt.*T.*Q, ... 84
               MoblFrac, ... 85
               MoblFrac.*T, ... 86
               MoblFrac.*T.^2, ... 87
               MoblFrac.*Q, ... 88
               MoblFrac.*Q.^2, ... 89
               MoblFrac.*T.*Q, ... 90
               TabltFrac, ... 91
               TabltFrac.*T, ... 92
               TabltFrac.*T.^2, ... 93
               TabltFrac.*Q, ... 94
               TabltFrac.*Q.^2, ... 95
               TabltFrac.*T.*Q] ; % 96
labels_SR4  = {'Const*';
               'T';'T^2'; ...
               'Q';'Q^2';'T*Q'; ...
               'PreScore*';'log(thetaE)';'log(thetaL)'; ...
               'D2';'D3'; 'D2*T'; 'D3*T'; 'D2*Q'; 'D3*Q'; 'D2*T^2'; 'D3*T^2'; 'D2*Q^2'; 'D3*Q^2'; 'D2*T*Q'; 'D3*T*Q'; ...
               'log(thetaE)*T'; 'log(thetaL)*T'; 'log(thetaE)*Q'; 'log(thetaL)*Q'; 'log(thetaE)*T^2'; 'log(thetaL)*T^2'; 'log(thetaE)*Q^2'; 'log(thetaL)*Q^2'; 'log(thetaE)*T*Q'; 'log(thetaL)*T*Q'; ...
               'pre*T';'pre*T^2';'pre*Q';'pre*Q^2';'pre*T*Q';'grade5';'grade5*T';'grade5*T^2';'grade5*Q';'grade5*Q^2';'grade5*T*Q';...
               'fem';'fem*T';'fem*T^2';'fem*Q';'fem*Q^2';'fem*T*Q';'blk';'blk*T';'blk*T^2';'blk*Q';'blk*Q^2';'blk*T*Q';'hsp';'hsp*T';'hsp*T^2';'hsp*Q';'hsp*Q^2';'hsp*T*Q';
               'SSEInc';'SSEInc*T';'SSEInc*T^2';'SSEInc*Q';'SSEInc*Q^2';'SSEInc*T*Q';'SSEIns';'SSEIns*T';'SSEIns*T^2';'SSEIns*Q';'SSEIns*Q^2';'SSEIns*T*Q';...
               '# HELPERS';'# HELPERS*T';'# HELPERS*T^2';'# HELPERS*Q';'# HELPERS*Q^2';'# HELPERS*T*Q';...
               'NoHmIntrnt';'NoHmIntrnt*T';'NoHmIntrnt*T^2';'NoHmIntrnt*Q';'NoHmIntrnt*Q^2';'NoHmIntrnt*T*Q';...
               'MoblFrac';'MoblFrac*T';'MoblFrac*T^2';'MoblFrac*Q';'MoblFrac*Q^2';'MoblFrac*T*Q';...
               'TabltFrac';'TabltFrac*T';'TabltFrac*T^2';'TabltFrac*Q';'TabltFrac*Q^2';'TabltFrac*T*Q'} ;
[alpha_SR4,se_SR4,tstat_SR4,pval_SR4,CI_SR4,out_SR4] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,[],[],0,[],[],[]); 

StDevEffectATE    = nan(length(alpha_SR4)+2,1) ;
Delta0 = alpha_SR4(1) + log(Y_EimputeShrunk)*alpha_SR4(8) + log(Y_LimputeShrunk)*alpha_SR4(9) ...
         + D2*alpha_SR4(10) + D3*alpha_SR4(11) + grade5*alpha_SR4(37) + fem*alpha_SR4(43) + blk*alpha_SR4(49) + hsp*alpha_SR4(55) ...
         + SSEInc*alpha_SR4(61) + SSEIns*alpha_SR4(67) + NHelpers*alpha_SR4(73) ...
         + NoHmIntrnt*alpha_SR4(79) + MoblFrac*alpha_SR4(85) + TabltFrac*alpha_SR4(91) ;
Delta0DescStats = [mean(Delta0) std(Delta0)] %#ok
StDevEffectATE(1) = std( Delta0 )/StDevDeltaScore ; %%%%This is a mean SD effect of Delta0

dDeltadT = (alpha_SR4(2) + D2*alpha_SR4(12) + D3*alpha_SR4(13) + log(Y_EimputeShrunk)*alpha_SR4(22) + log(Y_LimputeShrunk)*alpha_SR4(23) ...
              + prestdz*alpha_SR4(32) + grade5*alpha_SR4(38) + fem*alpha_SR4(44) + blk*alpha_SR4(50) + hsp*alpha_SR4(56) + SSEInc*alpha_SR4(62) + SSEIns*alpha_SR4(68) ...
              + NHelpers*alpha_SR4(74) + NoHmIntrnt*alpha_SR4(80) + MoblFrac*alpha_SR4(86) + TabltFrac*alpha_SR4(92) ) ...
            + 2*(alpha_SR4(3) + D2*alpha_SR4(16) + D3*alpha_SR4(17) + log(Y_EimputeShrunk)*alpha_SR4(26) + log(Y_LimputeShrunk)*alpha_SR4(27) ...
              + prestdz*alpha_SR4(33) + grade5*alpha_SR4(39) + fem*alpha_SR4(45) + blk*alpha_SR4(51) + hsp*alpha_SR4(57) + SSEInc*alpha_SR4(63) + SSEIns*alpha_SR4(69) ...
              + NHelpers*alpha_SR4(75) + NoHmIntrnt*alpha_SR4(81) + MoblFrac*alpha_SR4(87) + TabltFrac*alpha_SR4(93) ).*T ...
            + Q.*(alpha_SR4(6) + D2*alpha_SR4(20) + D3*alpha_SR4(21) + log(Y_EimputeShrunk)*alpha_SR4(30) + log(Y_LimputeShrunk)*alpha_SR4(31) ...
              + prestdz*alpha_SR4(36) + grade5*alpha_SR4(42) + fem*alpha_SR4(48) + blk*alpha_SR4(54) + hsp*alpha_SR4(60) + SSEInc*alpha_SR4(66) + SSEIns*alpha_SR4(72) ...
              + NHelpers*alpha_SR4(78) + NoHmIntrnt*alpha_SR4(84) + MoblFrac*alpha_SR4(90) + TabltFrac*alpha_SR4(96) ) ;
StDevEffectATE(2) = mean( BenchmarkChangeTstd*dDeltadT/StDevDeltaScore ) ;  %%%%This is mean pseudo-SD effect of T

dDeltadQ = (alpha_SR4(4) + D2*alpha_SR4(14) + D3*alpha_SR4(15) + log(Y_EimputeShrunk)*alpha_SR4(24) + log(Y_LimputeShrunk)*alpha_SR4(25) ...
              + prestdz*alpha_SR4(34) + grade5*alpha_SR4(40) + fem*alpha_SR4(46) + blk*alpha_SR4(52) + hsp*alpha_SR4(58) + SSEInc*alpha_SR4(64) + SSEIns*alpha_SR4(70) ...
              + NHelpers*alpha_SR4(76) + NoHmIntrnt*alpha_SR4(82) + MoblFrac*alpha_SR4(88) + TabltFrac*alpha_SR4(94) ) ...
            + 2*(alpha_SR4(5) + D2*alpha_SR4(18) + D3*alpha_SR4(19) + log(Y_EimputeShrunk)*alpha_SR4(28) + log(Y_LimputeShrunk)*alpha_SR4(29) ...
              + prestdz*alpha_SR4(35) + grade5*alpha_SR4(41) + fem*alpha_SR4(47) + blk*alpha_SR4(53) + hsp*alpha_SR4(59) + SSEInc*alpha_SR4(65) + SSEIns*alpha_SR4(71) ...
              + NHelpers*alpha_SR4(77) + NoHmIntrnt*alpha_SR4(83) + MoblFrac*alpha_SR4(89) + TabltFrac*alpha_SR4(95) ).*Q ...
            + T.*(alpha_SR4(6) + D2*alpha_SR4(20) + D3*alpha_SR4(21) + log(Y_EimputeShrunk)*alpha_SR4(30) + log(Y_LimputeShrunk)*alpha_SR4(31) ...
              + prestdz*alpha_SR4(36) + grade5*alpha_SR4(42) + fem*alpha_SR4(48) + blk*alpha_SR4(54) + hsp*alpha_SR4(60) + SSEInc*alpha_SR4(66) + SSEIns*alpha_SR4(72) ...
              + NHelpers*alpha_SR4(78) + NoHmIntrnt*alpha_SR4(84) + MoblFrac*alpha_SR4(90) + TabltFrac*alpha_SR4(96) ) ;
StDevEffectATE(4) = mean( BenchmarkChangeQstd*dDeltadQ/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD effect of Q

StDevEffectATE(7) = mean( ( alpha_SR4(7) + T*alpha_SR4(32) + (T.^2)*alpha_SR4(33) + Q*alpha_SR4(34) + (Q.^2)*alpha_SR4(35) + T.*Q*alpha_SR4(36) - Presig )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
StDevEffectATE(8) = mean( std(log(Y_EimputeShrunk))*( alpha_SR4(8) + T*alpha_SR4(22) + Q*alpha_SR4(24) + (T.^2)*alpha_SR4(26)+ (Q.^2)*alpha_SR4(28) + T.*Q*alpha_SR4(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE)
StDevEffectATE(9) = mean( std(log(Y_LimputeShrunk))*( alpha_SR4(9) + T*alpha_SR4(23) + Q*alpha_SR4(25) + (T.^2)*alpha_SR4(27)+ (Q.^2)*alpha_SR4(29) + T.*Q*alpha_SR4(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL)
StDevEffectATE(10) = mean( ( alpha_SR4(10) + T*alpha_SR4(12) + Q*alpha_SR4(14) + (T.^2)*alpha_SR4(16) + (Q.^2)*alpha_SR4(18) + T.*Q*alpha_SR4(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(11) = mean( ( alpha_SR4(11) + T*alpha_SR4(13) + Q*alpha_SR4(15) + (T.^2)*alpha_SR4(17) + (Q.^2)*alpha_SR4(19) + T.*Q*alpha_SR4(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3
StDevEffectATE(37) = mean( ( alpha_SR4(37) + T*alpha_SR4(38) + (T.^2)*alpha_SR4(39) + Q*alpha_SR4(40) + (Q.^2)*alpha_SR4(41) + T.*Q*alpha_SR4(42) )/StDevDeltaScore ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectATE(43) = mean( ( alpha_SR4(43) + T*alpha_SR4(44) + (T.^2)*alpha_SR4(45) + Q*alpha_SR4(46) + (Q.^2)*alpha_SR4(47) + T.*Q*alpha_SR4(48) )/StDevDeltaScore ) ;  %%%%This is a change from male to female
StDevEffectATE(49) = mean( ( alpha_SR4(49) + T*alpha_SR4(50) + (T.^2)*alpha_SR4(51) + Q*alpha_SR4(52) + (Q.^2)*alpha_SR4(53) + T.*Q*alpha_SR4(54) )/StDevDeltaScore ) ;  %%%%This is a change from wao to black
StDevEffectATE(55) = mean( ( alpha_SR4(55) + T*alpha_SR4(56) + (T.^2)*alpha_SR4(57) + Q*alpha_SR4(58) + (Q.^2)*alpha_SR4(59) + T.*Q*alpha_SR4(60) )/StDevDeltaScore ) ;  %%%%This is a change from wao to hispanic
StDevEffectATE(61) = mean( std(SSEInc)*( alpha_SR4(61) + T*alpha_SR4(62) + (T.^2)*alpha_SR4(63) + Q*alpha_SR4(64) + (Q.^2)*alpha_SR4(65) + T.*Q*alpha_SR4(66) )/StDevDeltaScore ) ;  %%%%This is a StDev change in SSEInc
StDevEffectATE(67) = mean( std(SSEIns)*( alpha_SR4(67) + T*alpha_SR4(68) + (T.^2)*alpha_SR4(69) + Q*alpha_SR4(70) + (Q.^2)*alpha_SR4(71) + T.*Q*alpha_SR4(72) )/StDevDeltaScore ) ;  %%%%This is a StDev change in SSEIns
StDevEffectATE(73) = mean( std(NHelpers)*( alpha_SR4(73) + T*alpha_SR4(74) + (T.^2)*alpha_SR4(75) + Q*alpha_SR4(76) + (Q.^2)*alpha_SR4(77) + T.*Q*alpha_SR4(78) )/StDevDeltaScore ) ;  %%%%This is a StDev change in NHelpers

OtherEffect    = nan(length(alpha_SR4)+2,1) ;
%%%%This is a mean total derivative of kappa with respect to Q, for a
%%%%change of BenchmarkChangeQstd, as a fraction of StDevDeltaScore
OtherEffect(4) = mean( ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd/StDevDeltaScore ) ;  %%%%This is a mean total derivative effect of Q (in original DeltaS units)
OtherEffect(7) = mean( ( T*alpha_SR4(32) + (T.^2)*alpha_SR4(33) + Q*alpha_SR4(34) + (Q.^2)*alpha_SR4(35) + T.*Q*alpha_SR4(36) )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score SLOPE EFFECTS ONLY (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
OtherEffect(8) = mean( std(log(Y_EimputeShrunk))*( T*alpha_SR4(22) + Q*alpha_SR4(24) + (T.^2)*alpha_SR4(26)+ (Q.^2)*alpha_SR4(28) + T.*Q*alpha_SR4(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE) SLOPE EFFECTS ONLY
OtherEffect(9) = mean( std(log(Y_LimputeShrunk))*( T*alpha_SR4(23) + Q*alpha_SR4(25) + (T.^2)*alpha_SR4(27)+ (Q.^2)*alpha_SR4(29) + T.*Q*alpha_SR4(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL) SLOPE EFFECTS ONLY
OtherEffect(10) = mean( ( T*alpha_SR4(12) + Q*alpha_SR4(14) + (T.^2)*alpha_SR4(16) + (Q.^2)*alpha_SR4(18) + T.*Q*alpha_SR4(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(11) = mean( ( T*alpha_SR4(13) + Q*alpha_SR4(15) + (T.^2)*alpha_SR4(17) + (Q.^2)*alpha_SR4(19) + T.*Q*alpha_SR4(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(12) = std( ( T*alpha_SR4(12) + Q*alpha_SR4(14) + (T.^2)*alpha_SR4(16) + (Q.^2)*alpha_SR4(18) + T.*Q*alpha_SR4(20) ) ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(13) = std( ( T*alpha_SR4(13) + Q*alpha_SR4(15) + (T.^2)*alpha_SR4(17) + (Q.^2)*alpha_SR4(19) + T.*Q*alpha_SR4(21) ) ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(37) = mean( ( T*alpha_SR4(38) + (T.^2)*alpha_SR4(39) + Q*alpha_SR4(40) + (Q.^2)*alpha_SR4(41) + T.*Q*alpha_SR4(42) )/StDevDeltaScore ) ;  %%%%This is a change from grade 6 to grade 5 SLOPE EFFECTS ONLY
disp('-----------------------------------------------------------------------------------')
disp('-------------------------  SHORT-RUN PROGRESS: MODEL 4  ---------------------------')
disp('-----------------------------------------------------------------------------------')
baselinetemp = sort( alpha_SR4(1)+prestdz*alpha_SR4(7)+log(Y_EimputeShrunk)*alpha_SR4(8)+log(Y_LimputeShrunk)*alpha_SR4(9) ...
                            +D2*alpha_SR4(10)+D3*alpha_SR4(11)+grade5*alpha_SR4(37)+fem*alpha_SR4(43)+blk*alpha_SR4(49)+hsp*alpha_SR4(55) ...
                            +SSEInc*alpha_SR4(61)+SSEIns*alpha_SR4(67)+NHelpers*alpha_SR4(73) ...
                            + NoHmIntrnt*alpha_SR4(79) + MoblFrac*alpha_SR4(85) + TabltFrac*alpha_SR4(91) - pre ) ;
BASELINE_SR4 = [mean(baselinetemp) median(baselinetemp) std(baselinetemp) baselinetemp(round(length(baselinetemp)*0.1)) baselinetemp(round(length(baselinetemp)*0.9))] %#ok
out_SR4.results.StDevEffectATE = StDevEffectATE;
out_SR4.results.OtherEffect = OtherEffect;
WhiteTest  = out_SR4.White %#ok
Method     = out_SR4.method %#ok
disp(out_SR4.results(:,:)) 
%%%%Test for joint significance of all intercept terms:           
R = zeros(16,length(ExpVars_SR4(1,:))) ; R(1,1)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; R(5,10)=1; R(6,11)=1; R(7,37)=1; ...
                                         R(8,43)=1; R(9,49)=1; R(10,55)=1; R(11,61)=1; R(12,67)=1; R(13,73)=1; R(14,79)=1; R(15,85)=1; R(16,91)=1; r = zeros(16,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
InterceptTest = temp.Wald %#ok
%%%%Test for joint significance of all NON-S^{pre} intercept terms:           
R = zeros(15,length(ExpVars_SR4(1,:))) ; R(1,1)=1; R(2,8)=1; R(3,9)=1; R(4,10)=1; R(5,11)=1; R(6,37)=1; ...
                                         R(7,43)=1; R(8,49)=1; R(9,55)=1; R(10,61)=1; R(11,67)=1; R(12,73)=1; R(13,79)=1; R(14,85)=1; R(15,91)=1; r = zeros(15,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
InterceptTest_NonSpre = temp.Wald %#ok
%%%%Test for joint significance of all T terms:           
R = zeros(48,length(ExpVars_SR4(1,:))) ; R(1,2)=1; R(2,3)=1; R(3,6)=1; R(4,12)=1; R(5,13)=1; R(6,16)=1; R(7,17)=1; R(8,20)=1; R(9,21)=1; 
                                         R(10,22)=1; R(11,23)=1; R(12,26)=1; R(13,27)=1; R(14,30)=1; R(15,31)=1; R(16,32)=1; R(17,33)=1; R(18,36)=1;
                                         R(19,38)=1; R(20,39)=1; R(21,42)=1; R(22,44)=1; R(23,45)=1; R(24,48)=1;
                                         R(25,50)=1; R(26,51)=1; R(27,54)=1; R(28,56)=1; R(29,57)=1; R(30,60)=1;
                                         R(31,62)=1; R(32,63)=1; R(33,66)=1; R(34,68)=1; R(35,69)=1; R(36,72)=1;
                                         R(37,74)=1; R(38,75)=1; R(39,78)=1; R(40,80)=1; R(41,81)=1; R(42,84)=1;  
                                         R(43,86)=1; R(44,87)=1; R(45,90)=1; R(46,92)=1; R(47,93)=1; R(48,96)=1;  r = zeros(48,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
TTest = temp.Wald %#ok
%%%%Test for joint significance of all Q terms:           
R = zeros(48,length(ExpVars_SR4(1,:))) ; R(1,4)=1; R(2,5)=1; R(3,6)=1; R(4,14)=1; R(5,15)=1; R(6,18)=1; R(7,19)=1; R(8,20)=1; R(9,21)=1;
                                         R(10,24)=1; R(11,25)=1; R(12,28)=1; R(13,29)=1; R(14,30)=1; R(15,31)=1; R(16,34)=1; R(17,35)=1; R(18,36)=1;
                                         R(19,40)=1; R(20,41)=1; R(21,42)=1; R(22,46)=1; R(23,47)=1; R(24,48)=1;
                                         R(25,52)=1; R(26,53)=1; R(27,54)=1; R(28,58)=1; R(29,59)=1; R(30,60)=1; 
                                         R(31,64)=1; R(32,65)=1; R(33,66)=1; R(34,70)=1; R(35,71)=1; R(36,72)=1; 
                                         R(37,76)=1; R(38,77)=1; R(39,78)=1; R(40,82)=1; R(41,83)=1; R(42,84)=1;  
                                         R(43,88)=1; R(44,89)=1; R(45,90)=1; R(46,94)=1; R(47,95)=1; R(48,96)=1;  r = zeros(48,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
QTest = temp.Wald %#ok   
%%%%Test for joint significance of all District 2 terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; R(1,10)=1; R(2,12)=1; R(3,14)=1; R(4,16)=1; R(5,18)=1; R(6,20)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
Dist2Test = temp.Wald %#ok
%%%%Test for joint significance of all District 2 slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; R(1,12)=1; R(2,14)=1; R(3,16)=1; R(4,18)=1; R(5,20)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
Dist2InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all District 3 terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; R(1,11)=1; R(2,13)=1; R(3,15)=1; R(4,17)=1; R(5,19)=1; R(6,21)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
Dist3Test = temp.Wald %#ok
%%%%Test for joint significance of all District 3 slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; R(1,13)=1; R(2,15)=1; R(3,17)=1; R(4,19)=1; R(5,21)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
Dist3InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District terms:           
R = zeros(12,length(ExpVars_SR4(1,:))) ; for ii=1:12; R(ii,9+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
DistTest = temp.Wald %#ok
%%%%Test for joint significance of all School District slope interactions:           
R = zeros(10,length(ExpVars_SR4(1,:))) ; for ii=1:10; R(ii,11+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
DistInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all School District Intercept Terms:           
R = zeros(2,length(ExpVars_SR4(1,:))) ; R(1,10) = 1 ; R(2,11) = 1 ; r = zeros(2,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
DistInterceptTest = temp.Wald %#ok       
%%%%Test for joint significance of all ThetaE terms: 
R = zeros(6,length(ExpVars_SR4(1,:))) ; R(1,8)=1 ; R(2,22)=1 ; R(3,24)=1 ; R(4,26)=1 ; R(5,28)=1 ; R(6,30)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
ThetaETest = temp.Wald %#ok           
%%%%Test for joint significance of all ThetaE interaction terms: 
R = zeros(5,length(ExpVars_SR4(1,:))) ; R(1,22)=1 ; R(2,24)=1 ; R(3,26)=1 ; R(4,28)=1 ; R(5,30)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
ThetaEInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL terms:          
R = zeros(6,length(ExpVars_SR4(1,:))) ; R(1,9)=1 ; R(2,23)=1 ; R(3,25)=1 ; R(4,27)=1 ; R(5,29)=1 ; R(6,31)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
ThetaLTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL interaction terms:          
R = zeros(5,length(ExpVars_SR4(1,:))) ; R(1,23)=1 ; R(2,25)=1 ; R(3,27)=1 ; R(4,29)=1 ; R(5,31)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
ThetaLInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all Student Type slope interactions:          
R = zeros(10,length(ExpVars_SR4(1,:))) ; for ii=1:10; R(ii,21+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
TypeInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; R(1,7)=1; for ii=1:5; R(ii+1,31+ii)=1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
PreScoreTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,31+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
PreScoreInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all age-cohort terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,36+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
Grade5Test = temp.Wald %#ok
%%%%Test for joint significance of all age-cohort slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,37+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
Grade5InteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all fem terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,42+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
FemTest = temp.Wald %#ok
%%%%Test for joint significance of all fem slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,43+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
FemInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all blk terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,48+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
BlkTest = temp.Wald %#ok
%%%%Test for joint significance of all blk slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,49+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
BlkInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all hsp terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,54+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
HspTest = temp.Wald %#ok
%%%%Test for joint significance of all hsp slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,55+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
HspInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEInc terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,60+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
SSEIncTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEInc slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,61+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
SSEIncInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEIns terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,66+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
SSEInsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEIns slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,67+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
SSEInsInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSE terms:           
R = zeros(12,length(ExpVars_SR4(1,:))) ; for ii=1:12; R(ii,60+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
SSETest = temp.Wald %#ok
%%%%Test for joint significance of all NHelper terms:           
R = zeros(6,length(ExpVars_SR4(1,:))) ; for ii=1:6; R(ii,72+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
NHelperTest = temp.Wald %#ok
%%%%Test for joint significance of all NHelper slope interactions:           
R = zeros(5,length(ExpVars_SR4(1,:))) ; for ii=1:5; R(ii,73+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
NHelperInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all Home Connectivity Controls:           
R = zeros(18,length(ExpVars_SR4(1,:))) ; for ii=1:18; R(ii,78+ii) = 1 ; end; r = zeros(18,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR4,'FGLSrobust',labels_SR4,R,r,0,[],[],[]); 
HomeConnectivityTest = temp.Wald %#ok
%%




dDeltadT = (alpha_SR4(2) + D2*alpha_SR4(12) + D3*alpha_SR4(13) + log(Y_EimputeShrunk)*alpha_SR4(22) + log(Y_LimputeShrunk)*alpha_SR4(23) ...
              + prestdz*alpha_SR4(32) + grade5*alpha_SR4(38) + fem*alpha_SR4(44) + blk*alpha_SR4(50) + hsp*alpha_SR4(56) + SSEInc*alpha_SR4(62) + SSEIns*alpha_SR4(68) ...
              + NHelpers*alpha_SR4(74) + NoHmIntrnt*alpha_SR4(80) + MoblFrac*alpha_SR4(86) + TabltFrac*alpha_SR4(92) ) ...
            + 2*(alpha_SR4(3) + D2*alpha_SR4(16) + D3*alpha_SR4(17) + log(Y_EimputeShrunk)*alpha_SR4(26) + log(Y_LimputeShrunk)*alpha_SR4(27) ...
              + prestdz*alpha_SR4(33) + grade5*alpha_SR4(39) + fem*alpha_SR4(45) + blk*alpha_SR4(51) + hsp*alpha_SR4(57) + SSEInc*alpha_SR4(63) + SSEIns*alpha_SR4(69) ...
              + NHelpers*alpha_SR4(75) + NoHmIntrnt*alpha_SR4(81) + MoblFrac*alpha_SR4(87) + TabltFrac*alpha_SR4(93) ).*T ...
            + Q.*(alpha_SR4(6) + D2*alpha_SR4(20) + D3*alpha_SR4(21) + log(Y_EimputeShrunk)*alpha_SR4(30) + log(Y_LimputeShrunk)*alpha_SR4(31) ...
              + prestdz*alpha_SR4(36) + grade5*alpha_SR4(42) + fem*alpha_SR4(48) + blk*alpha_SR4(54) + hsp*alpha_SR4(60) + SSEInc*alpha_SR4(66) + SSEIns*alpha_SR4(72) ...
              + NHelpers*alpha_SR4(78) + NoHmIntrnt*alpha_SR4(84) + MoblFrac*alpha_SR4(90) + TabltFrac*alpha_SR4(96) ) ;
dDeltadQ = (alpha_SR4(4) + D2*alpha_SR4(14) + D3*alpha_SR4(15) + log(Y_EimputeShrunk)*alpha_SR4(24) + log(Y_LimputeShrunk)*alpha_SR4(25) ...
              + prestdz*alpha_SR4(34) + grade5*alpha_SR4(40) + fem*alpha_SR4(46) + blk*alpha_SR4(52) + hsp*alpha_SR4(58) + SSEInc*alpha_SR4(64) + SSEIns*alpha_SR4(70) ...
              + NHelpers*alpha_SR4(76) + NoHmIntrnt*alpha_SR4(82) + MoblFrac*alpha_SR4(88) + TabltFrac*alpha_SR4(94) ) ...
            + 2*(alpha_SR4(5) + D2*alpha_SR4(18) + D3*alpha_SR4(19) + log(Y_EimputeShrunk)*alpha_SR4(28) + log(Y_LimputeShrunk)*alpha_SR4(29) ...
              + prestdz*alpha_SR4(35) + grade5*alpha_SR4(41) + fem*alpha_SR4(47) + blk*alpha_SR4(53) + hsp*alpha_SR4(59) + SSEInc*alpha_SR4(65) + SSEIns*alpha_SR4(71) ...
              + NHelpers*alpha_SR4(77) + NoHmIntrnt*alpha_SR4(83) + MoblFrac*alpha_SR4(89) + TabltFrac*alpha_SR4(95) ).*Q ...
            + T.*(alpha_SR4(6) + D2*alpha_SR4(20) + D3*alpha_SR4(21) + log(Y_EimputeShrunk)*alpha_SR4(30) + log(Y_LimputeShrunk)*alpha_SR4(31) ...
              + prestdz*alpha_SR4(36) + grade5*alpha_SR4(42) + fem*alpha_SR4(48) + blk*alpha_SR4(54) + hsp*alpha_SR4(60) + SSEInc*alpha_SR4(66) + SSEIns*alpha_SR4(72) ...
              + NHelpers*alpha_SR4(78) + NoHmIntrnt*alpha_SR4(84) + MoblFrac*alpha_SR4(90) + TabltFrac*alpha_SR4(96) ) ;



TotChng = ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd ; 
PartialChngQ = dDeltadQ*BenchmarkChangeQstd ; 
PartialChngT = dDeltadT*BenchmarkChangeTstd ; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%                 FIGURE 8:                 %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                               
figure
subplot(3,2,1) 
[Ftemp,Xtemp] = ecdf(PartialChngT) ;
stairs(Xtemp,Ftemp,'LineWidth',2)
xlabel('Partial Change $\frac{\partial\Delta\kappa_i}{\partial T_i}pSD(T)$ (in $36\!\times\!\kappa^{post}$ units)','Interpreter','latex','FontWeight','bold','FontSize',24)
ylabel('Empirical CDF','FontWeight','bold','FontSize',20)
xlim([-1.2,0.6]*StDevDeltaScore)
box on
grid on
subplot(3,2,3) 
[Ftemp,Xtemp] = ecdf(PartialChngQ) ;
stairs(Xtemp,Ftemp,'LineWidth',2)
xlabel('Partial Change $\frac{\partial\Delta\kappa_i}{\partial A_i}pSD(A)$ (in $36\!\times\!\kappa^{post}$ units)','Interpreter','latex','FontWeight','bold','FontSize',24)
ylabel('Empirical CDF','FontWeight','bold','FontSize',20)
xlim([-0.5,1.8]*StDevDeltaScore)
box on
grid on
grid on
subplot(3,2,5) 
[Ftemp,Xtemp] = ecdf(TotChng) ;
stairs(Xtemp,Ftemp,'LineWidth',2)
xlabel('Tot Change $(\frac{\partial\Delta\kappa_i}{\partial A_i}+\frac{\partial\Delta\kappa_i}{\partial T_i}\frac{\partial T_i}{\partial A_i})pSD(A)$ (in $36\!\times\!\kappa^{post}$ units)','Interpreter','latex','FontWeight','bold','FontSize',24)
ylabel('Empirical CDF','FontWeight','bold','FontSize',20)
xlim([-0.5,1.5]*StDevDeltaScore)
box on
grid on
subplot(3,2,[2,4,6]) 
[Ftemp,Xtemp] = ecdf(baselinetemp) ;
stairs(Xtemp,Ftemp,'LineWidth',2)
xlabel('Net Baseline Learning: $\Delta_{0i}-\kappa_{i}^{pre}$ ($36\!\times\!\kappa^{post}$ units)','Interpreter','latex','FontWeight','bold','FontSize',24)
ylabel('Empirical CDF','FontWeight','bold','FontSize',20)
xlim([-0.4,1.4]*StDevDeltaScore)
box on
grid on





%%
Q = AchievementData.Q ; Q = Q(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
T = AchievementData.T ; T = T(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS


pre = pretemp ;
diary off